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Abstract 

This proposal relates to the design, analysis and application of a novel numerical 
scheme for the solution of axisymmetric scattering problems. To this end, a procedure 
is introduced to iteratively evaluate the solution of the Lippmann-Schwinger integral 
equation in C'(Ai'log TV) operations, where is the number of the discretization points. 
The method achieves its efficiency through the use of the addition theorem and Fast 
Legendre Transforms (FLT). For globally smooth sound velocities/refractive indexes the 
method is spectrally accurate. More generally the order of convergence is tied to and in 
fact, limited by, the smoothness of the solution. 

1 Introduction 

A variety of numerical methods have been developed for scattering simulations. These include 
finite difference methods (FDM) (gS], gT], [HI), finite element methods (FEM) CpU]. fTT]. 
1^). integral equation methods (lEM) (|^, [H^, IS]) and methods based on fast Fourier 
transforms QOli g2]) and fast multipole expansions 

Each of these have advantages and shortcomings. For instance, FDM and FEM result in 
sparse matrices which can, therefore, be easily applied and inverted. Moreover, FDM can be 
implemented with relative ease. These methods however are typically tied to regular grids 
which prevents the attainment of high order convergence for arbitrarily shaped obstacles. 
In contrast, FEM, while more difficult to implement than its FDM counterparts can easily 
handle complex geometries, a feature that, in fact, largely justifies its popularity. A main 
disadvantage of the FEM approach (also shared by FDM), on the other hand, is related to 
the exterior nature of scattering problem. Indeed, this implies that the physical radiation 
conditions must be translated into conditions that allow for a finite computational domain. 
Although exact boundary conditions can be found they are non local (in space and time). 
For this reason, great effort has gone into the design of approximate local boundary conditions 
f|2U|. |5U| . |51) ) that minimize spurious reflections. To achieve this, however, the conditions 
must typically be imposed at a certain distance from the scatterer, leading to increased com- 
putational times and memory requirements. 



Integral equation methods, on the other hand, do not suffer from this problem, as they 
are based on integral formulas which implicitly account for radiation conditions through the 
use of outgoing Green's functions. These methods, however, typically lead to a linear sys- 
tem involving a full matrix and thus they arc not competitive unless a specialized strategy is 
used to accelerate matrix-vector products. Examples of accelerated lEM include those based 
on FFTs ((H21) ISD) and those that use fast multipole expansions (gS], EH)- Interestingly, 



1 



however, to date no implementation of these exists that can attain higher than first or second 
order of convergence. 

In this paper, we present a new integral equation approach to scattering simulations, 
especially designed to treat penetrable bodies of revolution. In the spirit of methods based 
on FFTs, our approach attains its efficiency through the use of fast Legendre transforms 
(FLTs). In contrast with classical methods, however, our scheme can be made to converge with 
arbitrarily large order. This important additional property is achieved using ideas inspired by 
work in jj] and [3] on two-dimensional scattering. As in 0] we use the addition theorem in the 
integral equation to derive a scheme whose convergence rate is tied to the global smoothness 
of the refractive index/sound velocity; for globally smooth characteristics, for instance, it 
converges spectrally. The contribution of this paper is to get the high order convergence of 
axisymmetric case in 3-D with two dimensional cost using FLTs. 



2 The Helmholtz and Lippmann-Schwinger equations 

As we mentioned, we are interested in the prediction of electromagnetic and acoustic scat- 
tering returns. Helmholtz equation arises naturally in acoustics and we shall concentrate on 
scattering problems for such models. Although our results will not be directly applicable to 
electromagnetics, they do suggest a strategy to follow in this case, as the corresponding inte- 
gral equations are similar in structure to those that arise in connection with the Helmholtz 
model. 

We will consider the case when the inhomogeneity is of compact support and the region 
under consideration is in K^. A bounded scatterer O is contained within a radius R from the 
origin. The refractive index n{x) varies arbitrarily in and n(x)=l outside the scatterer. 
Thus, setting m(x)=l — n{x)'^ we have m(x)=0 for x outside and a scattered field will 
be generated by an incident field satisfying the linearized equations in free space 

Au' + k^u' = in (1) 
The total field u = + u" will then satisfy 

Au + k'^n^u = in R3 (2) 

with appropriate conditions at infinity. The physically relevant condition is the Sommerfeld 
radiation condition which guarantees that the scattered field is outgoing, 

lim r ( ^ - iku'] = 0. (3) 

r^oo \ or J 

In fact, as we mentioned, this last condition constitutes one of the central challenges in 
the design of accurate numerical methods, especially if it must be approximated (see e.g., 
|47| . [24 ). Such approximations, however, can be avoided in an integral formulation of the 
problem that explicitly accounts for the outgoing character of the solution. Indeed, suppose 
u e C2(M3) is a solution of ©-© and let 

I gik\x-y\ 

^(x,y) = —- r, (4) 

^ ' 47r I a; - y I ' ^ ' 

denote the free-space outgoing Green's function. Let a; g IR'^ be an arbitrary point and choose 
an open ball B of radius r with exterior unit normal u containing the support of m and such 
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that X € B. Then, we arrive at the Lippmann-Schwinger equation which is equivalent to Q 
and lO, 

u{x) = - fc^ J^{x, y)u{y)m{y) dy, x e M.^ (5) 

B 

where B is any baU containing the support of m. 



3 Axisymmetric formulation of the Lippmann-Schwinger 
integral equation 

A convenient form of the Lippmann-Schwinger equation can be derived if we appeal to the 
addition theorem. The addition theorem (see states that 



I gik\x-y\ 

Att \ X — y \ 



ikY,hi^Hk\p\)Mk\p 



2n+ 1 



P„(cos6')F„(cos6i) 



(n — m)! 
(n + m)! 



P™(cos0)P™(cos6l) cos(to(0 - 0)) 



(6) 



where x — {psm9cosip,psiii9smip,pcos9), y — {p sin cos ip, p sin sin ip, p cos 0), and P™ 
are the associated Legendre functions (P^ — P„ are the Legendre polynomials). The series 
and its term by term first derivatives with respect to | a; | and | y \ are absolutely convergent 
on compact subsets of | a; |>| ?/ |. 

Therefore, we obtain, using the addition theorem and lemma, that 

r I gik\x-y\ 

<^{x,y)u{y)m{y)dy ^ —- -u[y)m{y)dy 



Att \ X — y 

P oo 

= tk lY.hlp{kpy)Mkp^) 



2n + l 
47r 



P„ (cos 6')P„ (cos 6*) 



(7) 



(n + m) ! 

rn—l ^ ' 

u ■ m ■ p^ sin dpdOdif, 



P^{cos0)P^{cos0) cos(m((^ - </.)) 



where 



P> = niax(p, p), p< = min(p, p). 
If the incoming ray and the obstacle are axisymmetric then, by uniqueness [H] , the solution 



M of (O and will be axisymmetric. Then simplifies to 

TT R 

• u{p, 0)m{p, 0)p^ sin dpd0. 
Finally, expanding u and in Legendre series we obtain 



R3 



^x,y)u{y)miy)dy^tkf^^^^^^ J j h\l\kpy)u{kp<)Pn{cos0)Pn{cos0) 
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n=0 n=0 



cosf j 

• e)m{p, 6)p^ sin 6* dpd6iP„ (cos d), (9) 

from which the axisymmetric formulation 

u^{p)^K{p)~ (^^+^^)'^' y j h^^\kpy)Ukp<)p^{cos9)u{p,e)m{p,e)p' sine dpde (lo) 



follows. 

4 Numerical implementation 

To implement pO|) numerically, we approximate u and m by a truncated Legendre series 



{p,cos8) = ^ u„(p)F„(cos6'), 



n=0 
A/ 



^(p, cos6') = ^ u„(p)P„(cos6'). 

n=0 

Then, if n < P, using the orthogonality properties of the Legendre polynomials we have 



u {p,cosO)Pn{cosd)m{p,cosd) sint 

{p,cose)Pn{cos0)^Tni{p)Pi{cos0) smOde, (11) 





2F 







1=0 

so that we m can be approximated by 

2F 



'(p, cos 61) = mi{p) Pi {cos 9). 

1=0 

The formulation l(Tn|) thus reduces to 



Un{p) = Kip) + '-K,,ip), (12) 



where 

R 



K„ip) = -{2n + l)fc3 / /i(i)(fcp>)j„(fcp<)/„(p)p2 (13) 







and 

TT 

/„(p) = Ju^{p,9)m^^{p,e)Pn{cos0)smede, (14) 



n = 0,1, . . . , F. To solve (|12f) . we shall appeal to an iterative method, e.g., GMRES, which 
reduces the overall problem to the fast evaluation of Kn{p) in H12|l . To achieve the latter 
we propose below an efficient strategy to calculate the angular integration in H14|) and the 
subsequent radial integration in (|13|l . 

4.1 Angular integration 

A natural approach to the evaluation of H14|l is to exploit the orthogonality of the Legendre 
polynomials by expanding the function m?^ {p,0) in Legendre series for fixed p. This, of 
course, demands that we calculate the Legendre transform {c„} 

■iF 

u^{p,e)m^^{p,9) = ^c„(p)P„(cos0), (15) 

ri=0 

which, in principle, requires 0{F^) operations (for fixed p). Fortunately, however, algorithms 
have been designed for the fast evaluation of the map 

{f{e,)}U ^ {^"}n=l, (16) 

where 

F 

/(^i) ^ c„F„(cos6'j), 



n=0 

and its inverse, for appropriate choices of the angles {9j}; see Appendix. These algorithms, 
which we shall refer to as Fast Legendre Transform (FLT) and Lnverse Fast Legendre Trans- 
form (IFLT), then suggest an implementation of the angular integration in the form 

{In{p)/rn}Lo = FLT3F{IFLnF{{un{p)}Lo) ' /FLTg^^ ({m„ (p) j^^o)) , (17) 

where 

1 

2 . 2 



Pn{x) I dx 



2n + l 



-1 



4.2 Radial integration 

As can be easily checked, the integrand in (|13|l has a corner-type singularity at p = p, which 
suggests that we write 



mill 



(a,R) R 



(27^^'l)fc3 =^n\ka) I Mkp)In{p)p^dp + Mka) I hi'\kp)In{p)p' dp 

min(a,R,) 
min(a,R,) R, 

Vnika) J jn{kp)In{p)p^ dp + jn{ka) j yn{kp)In{p)p^ dp 

min(a,R,) 

R 

+ jn{ka) j jn{kp)In{p)p^ dp. 





(18) 
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Since, hn\kp) blows up at the origin although jn{kp)hn\kp) behaves nicely. So, we define 
the modified Bessel functions jn{p), Vnip) as follows. 



jn{p) ■■= 



1- 3- 5 -...(271+1) 



Jnip) 



2P 



^n+l 



1 • 1 • 3 • 5 • . . . (2n - 1) 

Then, the equation (fT^ becomes 



Vn{p) 



l!(2n + 3) 2!(2n + 3)(2n + 5) 



l!(l-2n) 2!(l-2n)(3-2n) 



Kn{a) = i 



i(q,R) 



Vnika) 



,P\n+l 



R 



-jn{ka) 



min(a,R,) 



3n{kp)In{p)k^ pdp 



[-Tvn{kp)In{p)k'^pdp 

p 



R 



+ / 1.^3'^."5V'"(2n?i)^ ^"('^^'"(^^"'"^ 



(19) 



(20) 



and approximate each integrand in the right-hand side separately. To this end, we divide 
the integration domain in a number i 
^ ^ j ^ Ni on which we approximate 



the integration domain in a number Ni of equi-length interpolation intervals Uj = 



^^-1 1 

m=0 



(21) 



for p G Uj^ where 



p-(ui+u")/2 



H-^°)/2 

are Chebyshev polynomials in Uj. Then letting denote the Chebyshev points in Uj 
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we have 



p=l m=0 



v. 



/* 1 

m=0 ''^ 



ijn{kpi) 



Nd-1 



+ E 



m=0 



E E / (^)"y;(Mr:-"''(p)fc^pdp 



Pi 



-jn{kpi){kp] 



P=\ m = i-r 



(fcp)"(-2n- 
32-52...(2n+l)2 



(22) 



which can be readily evaluated if the moments 

pI+i 



/< 

Pi 

pi+i 



PY+'ukp)T:^^^hp)k^pdp, [ {^yY^{kp)T:^>hp)k'pdp, 



Pk 



{kpT 



l-3-5---(2n+l) 



Ukp)T::;^''^ {p)p' dp 



Pk 



[pQ = Up = Uj) are pre-computed and stored for 

< n < F, < m < Nd - 1, 1 < j < Ni, < k < Nd. 

4.3 Algorithm and operation count 

The prescriptions in §4.1 and §4.2 can be summarized in the following algorithm. 
Algorithm 

Input {uUpi)}n=o ' {™n(/'i)}n=o ■ Legeudrc transform coefficients ofu^{pl,0) andm{pl,0). 
Output {un{pi)}n^o ■' Legcndrc transform coefficients ofu{pl,9). 
Stages 
stage, 



• compute Bessel functions 

• compute ratios 



Vni^pl) , jnO^pi) 
' (4+1 )" 
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• compute moments aj^k.n.m, Pj,k,n,m, l3,k,n,m, 

pi+1 pk+1 . 

Ol3,k,n,m-= / {-^Y^^ in{kp)T^''^' {p)k'^ pdp , I3j^k,n,7yi ■= / { — )'^]M{kp)T"n''''' {p)k'^ P dp, 

J pk+1 l p 

Pi Pi 

7j,fe,n,™ := j i.-i.^...{2n + l) -'"^^' ^^'^ ^ 

Pi 

(e.g. with a Clenshaw — Curtis quadrature Iciyf ) 

• set Un{pi) = 

k stage, for k=l to the number of iterations required for GMRES to converge to a prescribed 
tolerance 

• compute angular integration {In{p'k)}n=o /"'^ given Un, 

{Wk)l^n}L, = FLTMlFLT,F{{un{pi)Y,^,) ■ IFLTM{^n{pi)}lL,)) 

• compute radial integration 

(a) compute {c^j^lTo^ 



/„(p)« clnT:^''^'{p) for peU, ■ 



(b) define p.j.k,n, Cj,k,n, Cj,fc/ 



Ni-l 

Pj,k,n ■= ^ ^ C?m'^j,k,n,m , Cj,k,n •= ^ ^ cinl^j,k,n,m , £,j,k,n ^ ^ C?ml j ,k ,n ,m] 

m— ?n— 772—0 
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(c) compute Kn{pl) 



p=l 1=0 



\n+l 



s{pZ-i) + Vn{kpZ.)pNuNd-l,n 



liPNd) = 3n{kpj^\)C,N,,Ni,n 



■q{pZa) + 3n{kp%''^_^)C,N,,Na-l,n 

,Nd-2,n 



i^n(pi):=*(s(pi)+g(pi))-fc^- 



3 • 5 • • • (2n - 1) 



(d) define 



As follows from the description above, at each iteration, the number of multiplications is 
given by 



Stage 


number of operations 


Angular integration 


0{F{\ogFfNANi) 


Radial integration 


(a) + (h) + (c) 


(a) 


OiNdilog Na)NiF) 


(b) 


0{NjN,F) 


(c) 


0{N,NdF) 



Therefore, the total number of operations per iteration is 

O [NiNdF{log'' F + log Nd + Nd + 1)] = 0{N log^ N) 
where N = NiN^F is the total number of degrees of freedom. 
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5 Numerical results 

In this section we present results that confirm the predicted performance of the algorithms 
described in § 4.3. 

Example 5.1 In our first example we consider the scattering off a homogeneous sphere of 
radius 1 and index n{x) — 2. In this case the problem is explicitly solvable, so that comparison 
with the exact solution is possible. Indeed, for a plane wave incident in the positive z-direction 
we have 

oo 

^ ^^kS{o,o^) ^ g^fcpcose = ^i"(2n+ l)j„(fcp)P„(cos6i). (23) 

n=0 

Equation \2S\) . together with the Helmholtz equation and the radiation condition (0) 
readily deliver 

^ iEZo{<^nJni2kp) - Z"(2n + l)j„(fcp)}^n(cOS0), p < 1 

where a„ and bn solve 

( -Jn(2fc) hn{k)\ ( a„\ _ _ (in{k)\ .X 

\~m2k) h'^{k)j[bj- Uwy^ 

Therefore, 



'Er=o{^P"j"(2M - l^ ,J%_,^ Jnikp)}Pn{cOs9), p < 1 

Xn=0 ( -1.1^.3^.'^"(2r-lV(2»+l) /^"j»(M + ibn^lM{kp))Pn{cOs9), p>l 



where solve 



a. 



jn{2k), -5n3n{k) - iyn(k) 

nj'„{2k) + 2kj'r!{2k), -Sninj'n{k)+j~n'ik)k)+iin + l)y;,{k) - iy;,'{k)kj \b 



'-^Jnikl 

:^inUk)+jJ{k)k)^ 



where 



-1- 12.32 •••(2n-l)2(2n+l) 

fl • 3 • • • • (2n- 1), n>l 
I 1 n^O 



In figure 1 we show the error E = Wu^xact — UapproxWoo between the approximate and exact 
solutions for different values of the interpolation orders Nd in Y21]) . The figure demonstrates 
that, indeed, the radial integrator incurs an error of 0{{jp^)^'^)- Numerical values corre- 
sponding to this figure are displayed in table 1 where we also exhibit timings and iteration 
counts. 
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-'Is 



1.2 1.4 1.6 1.8 
log,„(N|) 



2.2 2.4 



Figure 1: Error in radial integrator for Example 5.1 (sphere). 





time per iteration 


error 


log2 (error ratio) 


2-'' 


3 (sec) 


1.193720-03 




24 


8 (sec) 


6.51158e-07 


4.1963 


25 


15 (sec) 


3.91409C-08 


4.0563 


26 


30 (sec) 


2.41378C-09 


4.0193 


2^ 


61 (sec) 


1.50428e-10 


4.0042 



Table 1: Results corresponding to Example 5.1: sphere. 
Parameters: F = 2^ - 1, Nd = A, < p < 4. 

Example 5.2 To test the angular integration, our second example relates the sphere centered 
at (0,0, d) with radius r. Ui{p,6) is given as follows, 



yi _ gife(pcose-d) 



= 6-''='' ^i"(2n + l)jn{kp)Pn{cos( 



(24) 



n=0 



Since, 



Co i 

j Pm (cos 6) sin ede= j (t) dt 



V^(r^^^i2^P„i(cos(?o) 

fV(l-cos2 0o)g^i^^(cos0o), m>l 



where 



[{I- cos 6*0), 



m = 0. 



da;" 
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The exact Legendre series coefficients ofm{p, 6) is given as follows. ifO < {d—r) < p < {d+r), 
m{p,6) = 

(1 - ng)( £ V(l - cos^^o) ~ |j; P^(cos0o)Pr»(cos^) + i(l - cos0o)i'o(cos0) 

m=l ^ '' 

where (psin^Oi pcos^o) is a solution of 

+ = p^, + {z~ df = r^. 

So, we have the exact Legendre coefficients of u and m and an approximated solution is 
compared with the exact solution from example 1. Tables 2-5 shows the radial and angular 
convergences of a discontinuous scatterer. 





time per iteration 


error 


error ratio 


2^ 


1 (sec) 


0.0650681 




24 


2 (sec) 


0.0186698 


3.4852 


25 


3 (sec) 


0.00580343 


3.21703 


26 


7 (sec) 


0.00187212 


3.09993 


2^ 


13 (sec) 


0.00060713 


3.08355 


2« 


27 (sec) 


0.000201581 


3.01185 



Table 2: Results corresponding to Example 5.2: sphere centered at (0,0,2) with radius 1. 

Parameters: _F = 2^ - 1. N,, = 2. < p < 4. 





time per iteration 


error 


error ratio 


23 


2 (sec) 


0.0029432 




24 


3 (sec) 


0.00127203 


2.31378 


25 


7 (sec) 


0.000600751 


2.11741 


26 


14 (sec) 


0.000250397 


2.39919 


2^ 


27 (sec) 


8.41992e-05 


2.97387 


28 


54 (sec) 


1.89269e-05 


4.44865 



Table 3: Same parameters with table 2 except TV^ = 4. 





time per iteration 


error 


error ratio 


23 


4 (sec) 


0.000757108 




24 


7 (sec) 


0.000327485 


2.31189 


2-5 


14 (sec) 


0.000131013 


2.49964 


26 


27 (sec) 


3.06305e-05 


4.27719 


2^ 


53 (sec) 


7.49014e-06 


4.08945 



Table 4: Same parameters with table 2 except N4 = 8. 
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F 


time per iteration 


error 


log2 (error ratio) 


2^-1 


1 (sec) 


0.313367 




23-1 


1 (sec) 


0.0378805 


3.04833 


2^-1 


2 (sec) 


0.00503166 


2.91235 


2^-1 


5 (sec) 


0.00100595 


2.32247 


2^-1 


11 (sec) 


0.000180542 


2.47816 



Table 5: Angular convergence 
Parameters: iV, = 2^, = 2, 0<p< 4. 



Example 5.3 To test the convergence of the angular integration, our third example relates 
to the body of the revolution generated by an annulus. This hollowed sphere has the refractive 

index 

'(1 + VI cos 6* |P)i/2, l<p<2 
1, < />< 2 or p > 2. 



n{p, cos6) 



(25) 



Then, 



m{p,t) = 



^-\tf, 1<P<2 
0, 0<p<2orp>2. 



The Legendre coefficients ofm{p,t), 1 < p <2 is given by ^^-Qrn2n{p)P2n{t), 
where 



m2n{p) = -^^^ j P2n{t)\tfdt 

I 

' P2n{t)t'^ dt 



2n+ 1 



(4n + l)V7r2-^-ir(l 



(26) 



r(l-n+i/3)r(i/3 + 



r(i + /3) 



r(i + i/3)r(f + i/3) n;-^(i/3 + | + fc) 

In tables 6-9, we present the order of convergence for (3 =.2, 1.7, 2.2 and 3.2 respectively. 
The tables clearly show the correlation between smoothness of the refractive index n{x) and 
the order of convergence. 



F 


time per iteration 


error 


log2 (error ratio) 


2^-1 


1 (sec) 


0.000543237 


5.22994 


2^-1 


2 (sec) 


0.000123797 


2.13361 


2^-1 


5 (sec) 


2.75278e-05 


2.16901 


2^-1 


13 (sec) 


5.78823C-06 


2.2497 


28-1 


30 (sec) 


1.06212e-06 


2.44617 


29-1 


67 (sec) 


1.964336-07 


2.43484 
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Table 6: /? = .2, m e C•^ u e C^'^ 
Parameters: F = 2^° - 1, Nd = 2, A^i = 2^, < p < 4. 



F 


error 


log2 (error ratio) 


2^-1 


1.70623C-05 




2-^-1 


1.101666-06 


3.95306 


2*^-1 


8.65346e-08 


3.67025 


2'-l 


().(il41('-()!) 


.S.709G6 


28-1 


4.5934e-10 


3.84791 


29-1 


3.25004e-ll 


3.82104 



Table 7: /3 = 1.7, m e C^•^ u e C^'^ 
Same parameters with table 6. 



F 


error 


log2 (error ratio) 


24-1 


4.85002e-06 




2^-1 


1.87781C-07 


4.69086 


2^-1 


1.036716-08 


4.17896 


2^-1 


5.610646-10 


4.20771 


2S-1 


2.783076-11 


4.33342 


29-1 


1.409676-12 


4.30325 



Table 8: /3 = 2.2, m G C2■^ u € C^'^ 
Same parameters with table 6. 



F 


error 


log2 (error ratio) 


24-1 


3.435036-06 




2^-1 


5.821776-08 


5.88272 


2^-1 


1.326256-09 


5.45603 


2^-1 


3.581850-11 


5.21051 


2S-1 


8.99286-13 


5.31579 


2^-1 


2.39996-14 


5.22773 



Table 8: (3 = 3.2, m G C^ ^, u G C^'^ 
Same parameters with table 6. 
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6 Numerical analysis 

Actually, what we have solved is following, 

v{p,cose) = u'-^{p,cose) + ^{K^v){p,cose), (27) 



2 



where 



\p,cose) = ^<(p)P„(cos0) (28) 

ra=0 
F 

{K^v){p, COS 9) = Y.'^Kr,v){p)Pn{cose) (29) 



Therefore, 



While the exact solution satisfies, 



uM + ^iK^v){p), n<F 
0, n> F. 



u^{p,cos0) = M*'^(p,cos6i) + ^(i^^u)(p,cos6'), (30) 

Therefore, the error in the solution of the approximate integral equation H27() is given by 

\u{p,cos9) — {p,cos9)\ < {p,cosO) — {p,cos0)\ + \u^{p, cos0)\, (31) 

where 



T 

U 



{p, COS 9) = u{p, cos 9) — u^{p, cos 9) — m„(/3)P„(cos ( 



n>F 



By comparing (|?7jl and (pin)l . we get 



Therefore, 



F F " T^Fi F F\ I " T^F 1 
U — V — -K (U — V ) + -K U 



<B\\K^u^[ 

where, 



\{I--K^)-'\U<B. 



From the equations ifT^ . l(Tl|l . 



R ir 



{KnU^){p, COS 9) = -{2n+l)k^ J J h[l\kp>)in{kp<)p^{u^ {p,9)m'^^ {p,9)Pn{cos 9) sin 9 d9) dp, 



(32) 

Lemma 6.1 There is a constant C{R,k) such that 

R 

hi^Hkp>)jnikp<)p'dp\U. < -^M^ 

max(l, w^) 
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Proof. 



i(p,R) 



R 



I j h^n\kp>)jn{kp<)p'' dpiloo < \yn{kp)\ j b„(M|p' dp + \jn{kp)\ 



\y„{kp)\p^ dp 



min(/9,R) 



R 



+ 



\jn{kp)\ J \jn{kp)\p' 



dp. 



jn{p) 
Vnip) 

Therefore, 



hp' 



I '2 \ ■ 



l-3-5-...(2n+l) 
-l-l-3-5-...(2n-l) 



l!(2n + 3) 2!(2n + 3)(2n + 5) 



1 



\jn{ka)yn{kb)\ < 



l!(l-2n) 2!(l-2n)(3-2n) 



C{R,k) 



6"+H2n+ 1) 



Then, the result follows immediately. 



Therefore, 

\K^u^U<{2n+l)k\l\h^^\kp>)3^{kp^)\p'dp)Y, ^ ||mj(p)||oo||w«(p)||oo I \PmPsPn\dp 



R 



s+n 



< 



C{k,R) 



s>F ;=s— n 
1 



max{l 



-4tE E \\Mp)\\oo\\u,iP)\\oo \PmPsPn\dp 
s>Fl=s-n ■' 



where, 



m{p,cos6) = ^ mi{p) Pi {cos 9), 

m=0 

u^{p,cosO) = u s {p)Ps {cos 0), 

s>F 



(33) 



(34) 



Lemma 6.2 If f G C([— 1, 1]), the the Legendre coefficients of f, 



Cu = ^^^ I f{t)Pn{t)dx 



1 2n+l [ df(fy 
Itt 



n + l 2 



-{Pn-l{t)-tPn{t))dt. (35) 



Proof. Since, O = P„(cos^) satisfies, 



^(sin^^)=-n(n + l)(sin^)0, 
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Therefore, 



2n+ 1 



f{x)Pn{x)dx = / f{cos0)Pn{cose)smede 



' smb' — \dB 



n{n + 1) dO 



dd 



n{n + 1) 
1 



de 



1 



n{n +1) J de 



d If/ aw ■ Qd{Pnicos9)) 
— (/(cos&j) •smt' dO 



de 



d If, aw ■ adiPnicose)) 

■[J (cose)) ■ sme — de 



n{n + 1) J d9 



df{t) dt dPnit) dt 
n{n + l)J dt de dt M 



de 



dt 



dfjt) dPnjt) 
n{n + 1)7 dt dt 



(1 - t^)dt 



df{t) 



n+1 J dt 

-1 



{Pn-l{t) ~ tPn{t))dt 



For the last equaUty, we used the formula, — f^) — ?i(P„_i(i) — tPn{t)). 

Stieltjes' formula (01]) says, 

TT (2n + 1) a„ 



a2 



where 



(2 sin 6*) 2 
1 cos((n + |)0- ^) 
2n + 3 (2 sin 61) I 

1-3 cos((n + f)6>- ^) 

(27i + 3)(2n + 5) (2 sin 61) I 

_ (2n- 1)!! 
" (2n)!! 



(36) 



The series is convergent for ^ < < If the sum of the first m terms is taken as an 
approximation of P„(cos(?) for < 6* < tt, the error Pn,m{0) satisfies the inequality 



I I 4 1 1 
l^"'™' ^^(2n+l)a„ 



(2m - 1)!! 



(37) 



' {2n + 3)(2n + 5) . . . (2n + 2to + 1) (2 sin6l)'"+^ 



In other words, the error is less than twice the absolute value of the first term omitted if the 
cosine term is replaced by 1. 
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Lemma 6.3 If ^ is integrable, then the Legendre coefficients of f , 



\Cn\ = \^^ ff{x)Pn{x)\dx<-^ 



(38) 



Proof. Since P„(cos6') satisfies, 



4 11 
P„(cos0)-- — 
TT [Zn +1) an 



(2sin6l)^ 



4 

< -- 



1 1 



and, 



/ (cos 9) cos{ne) (sin B^- dO = f (cos 9) (sin B) 5 



sin(n6 



TT (2n+ l)(2n + 3) a„ (2 sin 

^ 1 J d{f {cos B){smB)i) 
n 



(i6i 



sin(n6l) dB 



Similarly, 



Therefore, 



/(cos 6*) sin(n6l) (sin 6*) ^ d6l = 0(i) 



'2n+ 1 



= 1 f{t)P.n{t)\dt 



< 



ficosB) 



4 1 1 



— 77; TT — COS ((n + - 

7r(2n+l)a„ 2 



1 , „ TT, (sin 61) 2 



4^ V2 



^6* 



7r(2n+l)(2n + 3)a„7 2i(sin6')5 




dB 



<o(4) 



where 



a„ w O(^) 



Corollary 6.4 If ^[pi- is integrable, then the Legendre coefficients of f , 



n 2 



(39) 



Proof. Combining 135|) and (|38|l . it follows directly. 



For example, if f{t) = g C'^+", < a < 1, then the Legendre coefficients decay 

like 0{ ^feVs ) by the corollary. But, numerically, it shows that f{t) decays as 0{ ^k+^s+o, )■ 
But, I will not attempt to prove the correlation between a and the convergence order here. 
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Lemma 6.5 Suppose {mi{p)}, \us{p)} from the equations iy4}) decay like C>(7Trm^); 0( ^k+i+a ) 
with a > respectively. If k > 1 then, 



If k — then, 



Proof. If A: > 1, 



\K U Hoc pk+i.b+a+min{l.b,k+a)'^' 



l^nU^lloo < V V ||mKp)||ooh.(p)||oo f\PmPsPn\dp 

^ ' > s>F l=s-n _j 

C{k,R) sr 1 1 }\p p p \ ,, 

- max{l,n) ^ ^ ;fc+2+a ,k+i+^ J \PrnPsPn\dp 



< 



s>F l=s — n 

C{k,R) ^1 11 



E 



max{l,n) {s — n)'^+i+" s'=+4-5+q ?naa:(l,n-^) 
^ C{k,R) 1 1 



Therefore, 



n 1 1 oo 
n=0 

F 



< 



C{k,R) 



-y 



Since 



n--5 (F + l-n)'=+^ 

71—1 ^ ^ 



1 1 



n=l ^ ' n=\ ^ 

^fc+Q ) ^ *^(^mm(fe+Q + .5,1.5) — ^mm(fe+c(, 1.5) 



Combining equations (|41|l and l|42() give our desired result. 
If fc = 0, 



s+n -'^ 



^ ' ' S>Fl = S-7l _^ 



< 



max{l, n) (F + 1 - n)-5+" F^-S+c 



(40) 



(41) 



(42) 



C{k,R) ^ 1 1 (43) 

C(A:,i?) 1 1 



19 



Therefore, 



tm n -r i 



C{k,R) 

F4-5+« max{l, n) (F + 1 - n)-5+« 



< 



1^1 \ £L I 

^ 1 ^ ]^ ]^ 



C(fc,i?) 

n=l ^ ' n=l ^ ' 



logF ^( log J- ) ^ ^( bg^^ 

' ' ^ pTnin{l,.5+a) ' ^p.5+a>' \ i:i4.5+r»-l-m,in,n ._5+r»W 



^4.5+a+min(l,.5+a) - 

(44) 



Lemma 6.6 If {mi{p)}, {us{p)} decay like O(pE^), 0( ^^+2+0 ) respectively where k = 0, 1. 
Then, 

i0{-^) ifk=l, 



F„,T\\ 



\\K^u 



\0{t^) ifk = 0. 



Proof. If A; = 1, 

' s>Fl=s-n Zi 



Similarly, 



Therefore, 



If A; = 0, then 



s>F 



{s - n)i+« s3-5+a maa;(l,n-5) 



< 



C{k,R) 



max{l, n-5) (F + 1 - n)« F^.s+c 



n-5 (F + 1 - n)« 



< 



E 



= 0( 



J7'-.5+a 



1 



^min( — .5+a,.5) 



1 



(45) 



(46) 



(47) 



^ ' s>Fl=s-n _i 
<C{k,R)y ^ 

C{k,R) 



PmPsPn\dp 



1 



< 



(s - n)" s2.5+a TOax(l, n-5) 

1 1 



maa;(l, n-^) (F + 1 - n)" 



(48) 
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Therefore, 



F„,T\ 



\K''u 



0{ 



pl+2c 



(49) 



Corollary 6.7 If{mi{p)}, {us{p)} decay like 0{:^^^), 0( ^^+2+0 ) respectively. Then, 

0{ pk + 2.5 + a + min(1.5,k~2 + a) ) if k > 3 



Sf := \\UF ~ vf\ 

Mf^) 

From (|31() . the total error is given as 

\\u{p,cos9) — cos6')||oo < \\u^ {p,co89) — t;^(p, cos6')| 

where 



ifk = 2, 
ifk^O. 



\\u^{p,C08t 



(50) 



{p,cos9) — u{p, cos 9) — {p, cos 9) — ^ w„(p)P„(cos( 



n>F 

So far, we have error bound for ||wf ^ Wf^||oo and the error bound for ||u^(p, cos^)||oo is easily 
given as 

CTF II Unip)Pn{cOs9)\\oo = Q( ^fc+i+„ ) 

n>_F 

Where u„(p) ~ (^( ^fc,,,^^^ ) Therefore, the total error is dominated by 

If m{p,9) is given as then, m„(p) ~ ^^^^^ and m„(/o) ~ 0( „2.5+)!) )- Therefore, the 
total error is given as 



1 



||u(p, cos6') - -y^(p, cos6')||oo ~ 0{ 

. But, actually, our tables 6-9 show ||m — w^Hoo ~ Q( j^2+g )■ 

The reason is that <j_F — \\ X],i>_f Wn(p)-Pn(cos 0)||oo converges faster than tpF ■= \\ J2n>F l^n(p)lll 
and the latter is 0{ j^ii+^j )■ For example, when f3 = 2.2, 



F 


total error 




l0g2('^F) 


ifF 


log2(<^F) 


24-1 


4.85002e-06 


4.82746e-06 




2.4488e-05 




2^-1 


1.87781e-07 


1.87812e-07 


4.6839 


1.72411e-06 


3.82815 


2^-1 


1.03671e-08 


1.03674e-08 


4.17917 


1.28051e-07 


3.75106 


2^-1 


5.61064e-10 


5.61065e-10 


4.20774 


9.72137e-09 


3.71942 


28-1 


2.78307e-ll 


2.78301e-ll 


4.33345 


7.44753e-10 


3.70633 


29-1 


1.40967e-12 


1.40934e-12 


4.30355 


5.29085e-ll 


3.81519 



Table 10: f3 = 2.2, m e ^ ueC 



Parameters: F 



1, Nd = 2, N, = 2\ < p < 4. 
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A Appendix 



A.l The Fast Legendre Transform and its Inverse 
A. 1.1 General theory of orthogonal polynomials 

Lemma A.l (three-term recurrence) Let {pm} be an orthogonal polynomial sequence for a 
nonnegative integrable weight function. Then {pm} satisfies a three-term recurrence relation 

Pk+i{x) = (Akx + Bk)pk{x) + CkPk~i{x), (51) 

where Ak, Bk, Ck are real numbers with Ak 7^ and 7^ 0. 

Proof. See, e.g., [0, Theorem 4.1] ■ 

Next define the associated polynomials Qi^rm Ri,m for the orthogonal polynomial se- 
quence {pm} by the following recurrences on m ^ 

Qum{x) ^{Ai+m-lX + Bi+m-l)Ql,m-l{x) + Cl+rn-lQl,m-2{x) , 

Qi,o{x) = 1, Qli = Aix + Bi, ^^^^ 

Rl,m{x) =(Ai+m_ia; + Bi+m-l)Rl,m-l{x) + Cl+m-lRl,m-2{x) , 

i?,,o(x) -0, Ri,i^Ci. 

We have 

Lemma A. 2 (generalized three-term recurrence). The associated polynomials satisfy degQi = 
m, deg Ri^m < m — 1, and for I > 1 and m > 1, 

Pl+rn = Ql,m ' Pi + Rurn ' Pl-1- (53) 

Proof. Equation H53|) follows by induction on m with the case m = 1 being the original 
three-term recurrence (|51f) . ■ 



A. 1.2 The Legendre Transform 

Assume that f{x) is a polynomial of degree less than and let 

N-l 



where 



and 



Then, 



/(^) = X! c„T„(a;). 



(-l)"£n , (2j + l)nn 
c» 2. cos — (54) 



1, ifn = 0, 

2, if n > 0. 



N-l 

I f{x)dx= 1 /(cos e) smOdO^^, 

i n=0 



AT-l 
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where 

Then with (|54|l we obtain, 

/f^ W (2z + l)7r (2z + l)7r 

/ = }^ ^ a„ /(cos ) cos(n z) 

n=0 i=0 

i=0 n=0 

N-1 

/(cos — — — )w, = 2^ fix, )w, , 

i=0 1=0 

where 

J, (-l)"e„ (21 + 1)71 

= -^^anCOs(n ). 

In particular, if the degree of f{x) is less than or equal to N, then 

J. 2^-1 

Pmix)fix)dx^ Pmix^'')f{xDwr, m G 0, 1, 2, , TV _ 1. 
Definition A. 3 The Legendre coefficients {cm} of f{x) are defined as 

2N-1 

c^ = ±J2 P™(xr)/(^rK'^, meO,l,2,,iV-l, 

1=0 

1 

where {Pm} ifc the Legendre polynomials and Tm — / | Pm{x) \ dx = 2m+i ■ define the 

-1 

Legendre Transform as the map 

{fi^i )}i=0 ^ {cm}„i=o- 

A. 1.3 The Fast Legendre Transform (FLT): Healy-Driscoll algorithm [.14] 

In [33], Inda, Bisseling and Maslen reviewed and implemented an original idea of Healy- 
Driscoll 55 to compute 

N-l 

J2 P,n{xf)f{x^), m e {0, 1, 2, , AT - 1} 

i=0 

in 0{N log^ N) operations in the case where {pm} is any sequence of orthogonal polynomials, 
xf are Chebyshev points and N is power of 2. 

Here we recall the ideas behind the scheme. First we define the truncation operator 7^ 

by 

n-l 

%f=Y.bkTk, (56) 

fe=0 



TT 

an— J cos{n9) sinO d9. 
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where 

oo 

and the truncated polynomial by 

Zf" = rK{f -pi), (57) 

where {pi} are orthogonal polynomials. We also let 5^ denote the Lagrange interpolation 
operator Sn that is, Snf{x) is the polynomial of degree less than n which agrees with f{x) at 
the Chebyshev points ccq , . . . x^_i. 

If the degree of f{x) is less than m, then for n < to we have 

^n/(x) = ^6fcTfe, (58) 

fc=0 

where 

{&fc}r=o'=FCT({/(4')};ro'). 
Therefore, for an arbitrary polynomial ,f{x) and n < m, 

n-l 

r„(5„/(a:)) = ^&fcrfc, (59) 

fc=0 

where 

{&fe}r=o'=FCT({/(a;7)}™-o'). 

When the degree of the polynomial f{x)pi{x) is less than 2N, due to the Gaussian quadrature 
(??) we have 

Zl^T,if-p,)^^jl^^dx^^j:p,{x^)f{xr). (60) 

Therefore a fast Legendre transform will follow from a scheme that computes Zf fast, 
which is precisely what the Healy-Driscoll algorithm is designed to do. To introduce this 
procedure, we first note that, from (|53|l . 

/ • Pl+K = Ql.K ■ if ■ Pi) + Rl,K ■ if ■ Pl-l), (61) 

and therefore 

TkU ■ PI+k) = TKiQLK ■ if ■ Pi) + Rl.K ■ if ■ Pi-i)). (62) 
On the other hand, from (??), 

Qi,KT2K+i S spanj{{T/f+i, . . . , T^K+i}, 

Ri,KT2K+i G spanjj{TR-+i+i, . . . ,T^K+i-i}- 

Therefore, 

rK{,Ql,K ■ if ■ Pi) + RlK ■ if ■ Pl-l)) - TKiQLK ■ TkU ■ Pi) + Rl,K ' ^Kif ■ Pl-l))- (63) 
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The degree of Qi^k • T2K if ■ Pi) is less than 3K and by applymg the Lagrange interpolation 
operator, we obtain 



■ 'hnif -pi) - S2k{Qi,k ■ TikU ■ Pi)) = T2K ■ w{x), 



where w{x) is polynomial of degree less than K and we have used the fact that the Chebyshev 
points x'j^ are the zeros of T2K ■ 



Again, from (??), 
Therefore, 



T2K ■ w{x) e spanjj{T/f , . . . , T^k}- 

Tk{T2k-w{x)) = Q. (64) 



(65) 



Finally, from 

Z^+K = Tk{Qi,k ■ if ■ Pi) + Ri.K ■ if ■ Pl-l)) 

= TxiQl^K ■ T2Kif -Pl) + Rl.K ■ TlKif -Pl-l)) 
= TKiS2KiQl,K ■ T2Kif ■ Pl) + RlK ' r2Kif ' Pl-l))) 
= TR'(52if (^/ • Qi,k) + S2KiZi_i ■ Ri^k))- 

With a similar argument, we can also show that 

Zl+K-l = TKiS2KiZi^ ■ Ql,K-l) + S2KiZ^^i ■ Rl^K-l))- (66) 

The Healy-DriscoU algorithm is based on formulas (|65|l and 16611 . From (|59|) . Z;^^ -^, 
Z^j^ can be computed from Zj^^ and Zf^ using the FCT in OiKlogK) operations. At 
stage 0, the algorithm computes (Zq, Z-[) using FCT. At stage 1, we get iZ^'^l, Z^Jj^^^) 
from iZ^ , Zi) using (|()^ . (|?^ . At stage 2, we compute (Z^|^, Z^j'^_^^) from iZQ^^, Z^^^) 
and (Z^^^, ^^4^_|_i) from iZ^j^, figure 2. The total number of stages 

is log2(A'') and each stage costs 0(7Vlog7V) operations, for an overall cost of 0(7Vlog^ A^) 
operations. 
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1 3 5 7 9 11 13 15 

Figure 2: Computation of Zj^ with = 16 

A. 1.4 The Inverse Fast Legendre Transform (IFLT) 

Clearly the IFLT can be written as a matrix in the form 



M 



N 



As shown in |21j . 



\Pn-i{xq) pn-i{xi) 



PN~l{x^) \ 
PN-liXi) 

PN-lix%_i)J 



Po{x%^i) \ 

Pli^N^l) 
PN-l{x%_i)J 



can be applied to a vector in 0(A^(log A^)^) operations. Indeed can be decomposed into 
log2 N factors, 

M% = Aiog^ nAios^ n~i--- A2A1 (67) 

(68) 



and A\ can itself be applied in 0(iV(logiV)) operations. Therefore we have that 

IFLTm = Mn = A\A\ - ■ ■ Afog^ w-i Aog2 AT 
can be applied in 0{N\o^ N) operations; for details, see |21| . 
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Lemma A. 4 IFLT algorithm. 

Input f = (/o, . . . , /]v_i); Legendre transform coefficients and N is a power of 2. 

Output f = (/o, . . . , /at-i); Inverse Legendre transform values. 

Stages 

stage, 

a=[/o,0,/i,0,...,/;v-i,0] 

where ai = [fi, 0] and a = [ao, • • • , o-n-i] 
k stage, for k=l iolog2 N — 1 

K = 2'' 

(a) calculate c e A^i*4jv, 

c = [cq, . . . ,C2N/K-i], Ci = cheby ^Kiio-i^oi^K]), Ci e Mu2K 
where oi*k '■ zero matrix of size 1 * K 



(b) calculate g e Mu2n, 
(A^ O 



O A^' 



\o o 



o 



( 4 \ 



NI{2K)J 



( 4 \ 



9\ 



Qi e Mu2K 



\'^2N/K-1/ 



\9n/k-i/ 



where 



K ( I2K*2K 02K*2K R, 



02K*2K ^2K*2K Qi 



Rf-^ = diag {Rf-\xl''), iif "'(arif e M2k*2k 
Rf= diag{Rf{xl''),...,Rf{xl^_,)) 
Qf-i = dtag {Qf-\xl^), Qf-\xl^_,)) 
Qf= diag{Qf{xl^ ),..., Qf[xll_^)) 

; = * (i - 1) + 1 

(c) calculate a € M.\*2n, 



a= [ao,...,ajv/i<r-i], a, = cheby 2K{9i) 



• Given a = [ao, ai], compute e e A^i*2JV 

e = [ icheby]^{ao), ichebyf^{ax) 

• f = (/o, . . . , /jv-i) «5 computed as 



{diag{Po{xf)) diag{Pi{xf))) (e*) = (f*) 
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A. 1.5 Numerical examples 



In table 5 and figure 3 we present results representative of the performance of the algorithms 
described in A. 1.3 and A. 1.4. To assess the stability and accuracy of the algorithms we shall 
compare their outcome to direct calculations of the transform and its inverse. To this end, 
we let 

i?LT({/(.r)}S,-^) - Ml-Jo (69) 
denote the forward transform, and 

IDLTiMZzl) = {/(o^DlS"' (70) 

the inverse, computed directly, without acceleration and using the three-term recurrence ()51|l 
to evaluate Legendre polynomials. We further introduce the notation DLT'^^ and IDLT'^^ 
to denote the same transforms but with the three-term recurrence ()51|l evaluated in quadruple 
precision. In table 5 we compare the errors 

\\DLTi,{IDLT^''{v'^'^)) - v'^'^U, 

\\IDLTn{v^^') - /i?L^^^(^;f ')|U 
labeled "DLT" and "IDLT" respectively, with that incurred by the fast transforms, namely 

WFLTilDLT^^'iv^/')) -v^/^Woo ("FLT") 

\\IFLT{v^/^) - IDLT^''iv^/')\\^ ("IFLT") 

where is a vector of size m with unit nth component and zeros everywhere else. In ad- 
dition, the table includes results obtained on application of the fast transforms when the 
recurrence (IS^ are precomputed in quadruple precision ("FLT-QP" and "IFLT-QP"). We 
see that this latter strategy provides errors that are comparable to double precision DLT and 
IDLT calculations. Finally, in figure 3, we display timings for the FLT algorithm showing 
that, indeed, the operation count is proportional to iVlog^ N. 



N 


DLT 


FLT 


FLT-QP 


IDLT 


IFLT 


IFLT-QP 


256 


1.98476e-14 


2.8107e-13 


1.17584e-14 


2.52909e-13 


5.27689e-13 


1.99396e-13 


1024 


8.0331e-14 


3.49153e-12 


3.21525e-14 


4.18487e-12 


2.3919e-10 


4.22617e-12 


4096 


1.49405e-13 


1.02863e-10 


1.51153e-13 


8.09853e-ll 


1.45763e-08 


3.95126e-ll 


16384 


9.05882e-13 


1.27641e-09 


5.2725e-13 


1.2968e-09 


4.37189e-07 


5.39e-10 



Table 5: Examples of errors in the FLT and IFLT algorithms, QP indicates that the (pre-) 
computation of the recurrences H52|) are carried out in quadruple precision. 
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Figure 3: Timings for the evaluation of FLT({/(.Tf'' )}|iV ) 
and lFLT{v^^'^) for different values of N. 

A. 2 Computation of modified Bessel functions 

The spherical Bessel functions j^{p),y,^{p),h^\p),hl^\p) are defined in terms of 
Bessel functions by the relations, 



hi'^^^'\p) 



where 



h(^^^^^\p)=Mp)±iy.{p)- 

A standard calculation of Bessel functions evaluates these from the recurrence 

Zu+i{p) = — Z^{p) - Z^_i{p), 
P 

Our modified Bessel functions are as follows. 



^ 1 ■ 3 ■ 5 ■ ■ ■ ■ (2n + 1) . 

3nKP) = — Jn[p) 



1 



+ 



Cip'r 



l!(2n + 3) 2!(2n + 3)(2n + 5) 



+ ... 



Vnip) 



r2/n(p) 



-1-1-3-5 - ...(2n- 1) 

Therefore the new recurrences of these are 



1 - 



hp' 



+ 



iip'r 



l!(l-2n) 2!(l-2n)(3-2n) 



jn-l(p) =jnip) 



(2n + l)(2n + 3) 



jn+i{p) 
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where, 



p2 ^ 

yn+l{p) = Vnip) - -77, —:Vn-l{p) (79) 

(2n — l)(2n + Ij 



yo{p) = -P- \ I^Yi{p) 



2p 



mip)^-p'-^Yp^i{p) 

The reason that we need to compute j„ (p) with downward recurrence is that (|78|l has 2 hnearly 
independent solutions jnip) and j/n(p) = "'^^^ '^l^""'"^'' yn{p)- Since j/n(p) is an exponentiaUy 
growing solution as n increases, the upward recurrence relation is numerically unstable since 
the round-error in the jo{p) and ji(p) is rapidly amplified by the recurrence. On the other 
hand, the downward recurrence is stable since the round-off error in the starting values is 
quickly damped by the recurrence. Computing jn{p) and jn~i{p) can be done directly from 
the f7^. 
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